function FPU_cal()
    N = 5;
    K = 1;
    omega0 = K/1;
    alpha = 0.0;
    t0 = 1e2; %单位为us
    
    %体系状态初始化
    ini_state = cell(1,N);
    for i = 1:N
        ini_state{i} = [0,0];
    end
    ini_state{1} = [1,1];
    en0 = tot_energy(N, ini_state, K, alpha); %计算体系的初始能量
    
    %设置所得数据的参数
    tau = t0 * 1e-6; %Fourier级数系数的采样间隔
    %sam_len = 1e3; %Fourier级数系数的样本长度
    step = 1e3; %演化的步长，同时也是Fourier级数系数的样本长度
    list_len = 1e4; %总演化数
    t_min = step * t0 * 1e-6;
    t_max = step * list_len * t0 * 1e-6;
    time_list = linspace(t_min, t_max ,list_len);
    
    %设置输出参数
    %error_list = zeros(1,list_len);
    energycell = cell(1,N); %记录各模式的能量随时间的变化
    xcell = cell(1,N); %记录各质点位移随时间的变化
    xtcell = cell(1,N); %记录各质点速度随时间的变化
    %omegacell = cell(1,N); %记录简正频率（1到N 对应 频率从小到大）随时间的变化
    for i = 1:N
        energycell{i} = zeros(1,list_len);
        xcell{i} = zeros(1,list_len);
        xtcell{i} = zeros(1,list_len);
        omegacell{i} = zeros(1,list_len);
    end
    
    %计算alpha=0时的本征模式和本征频率，以确定寻找位移本征模式的傅里叶级数系数峰的最小峰值
    [evector, omega] = eigen_cal(N);
    omega_max0 = max(max(omega));
    
    %用Fourier级数展开计算本征模式和本征频率
    
    
    state = ini_state;
    
    for j = 1:list_len
        
        statej = v_verlet(N, state, t0, step, K, alpha);
        state = statej;
        
        en = tot_energy(N, statej, K, alpha);
        
        %error_list(j) = abs(en-en0)/en0; 
        error_list(j) = en;
        
        for k = 1:N
            xcell{k}(j) = statej{k}(1);
            xtcell{k}(j) = statej{k}(2);
        end
        
        xlist = zeros(1,N);
        xtlist = zeros(1,N);
        for k = 1:N
            xlist(k) = statej{k}(1);
            xtlist(k) = statej{k}(2);
        end
        
        for k=1:N
            qcell{k}(j) = xlist * evector(:,k);
            qtcell{k}(j) = xtlist * evector(:,k);
        end
        
        for k=1:N
            energyj = 0.5 * (qtcell{k}(j))^2 + 0.5 * ((omega(k,k)/omega0)^2) * (qcell{k}(j))^2;
            energycell{k}(j) = energyj;
        end
        
    end
    
    mode_en_error_list = zeros(1,list_len);
    for q = 1:list_len
        en_mode = 0;
        for r = 1:N
            en_mode = en_mode + energycell{r}(q);
        end
        %mode_en_error_list(q) = abs(en_mode - en0)/en0;
        mode_en_error_list(q) = en_mode;
    end
    
    %将数据存为txt文件
    fname1 = sprintf('N=%d_K=%.1f_alpha=%.3f_t0=%dus_step=%d_en&x&xt.txt',N, K, alpha, t0, step);
    fid1 = fopen(['.\',fname1],'w');
    for m = 1:list_len
        fprintf(fid1, '%d\t',time_list(m));
    end
    fprintf(fid1,'\n');
    for m = 1:list_len
        fprintf(fid1, '%d\t',error_list(m));
    end
    fprintf(fid1,'\n');
    for m=1:N
        for n=1:list_len
            fprintf(fid1, '%d\t',xcell{m}(n));
        end
        fprintf(fid1, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid1, '%d\t',xtcell{m}(n));
        end
        fprintf(fid1, '\n');
    end
    fclose(fid1);
    %%%%%%%
    fname2 = sprintf('N=%d_K=%.1f_alpha=%.3f_t0=%dus_step=%d_modes-en&q&qt.txt',N, K, alpha, t0, step);
    fid2 = fopen(['.\',fname2],'w');
    for m = 1:list_len
        fprintf(fid2, '%d\t',time_list(m));
    end
    fprintf(fid2,'\n');
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',qcell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',qtcell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    for m=1:N
        for n=1:list_len
            fprintf(fid2, '%d\t',energycell{m}(n));
        end
        fprintf(fid2, '\n');
    end
    fclose(fid2);
    
    
    %画图
    figure(1)
    scatter(time_list, error_list,'DisplayName','x');
    hold on 
    scatter(time_list, mode_en_error_list,'DisplayName','q');
    xlabel('time(s)');
    ylabel('energy');
    legend('show');
    
    figure(2)
    for p=1:N
        dotsname = sprintf('particle%d',p);
        plot(time_list, xcell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('x(a.u.)');
    legend('show');
    
    figure(3)
    for p=1:N
        dotsname = sprintf('particle%d',p);
        plot(time_list, xtcell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('v(a.u.)');
    legend('show');
    
    figure(4)
    for p=1:N
        dotsname = sprintf('mode%d',p);
        plot(time_list, qcell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('q(a.u.)');
    legend('show');
    
    figure(5)
    for p=1:N
        dotsname = sprintf('mode%d',p);
        plot(time_list, energycell{p},'DisplayName',dotsname);
        hold on
    end
    xlabel('time(s)');
    ylabel('energy(a.u.)');
    legend('show');

end